% 参数设置
global G m
G = 1; % 引力常数
m = 10000; % 中心天体质量
m_o = 1; % 运动质点质量
E = [-5, 0, 5]; % 总能量，分别对应E<0, E=0, E>0
L = [0.5, 1, 2]; % 角动量
tspan = [0, 100]; % 时间范围
r0 = 100; % 初始位置
v0 = sqrt(2*E + 2*G*m*m_o/r0); % 初始速度，调整以匹配不同的E值
line(0,0,'color','r','marker','*')
for i=1:3
    [t,u] =ode45 (@fffunction, 0:0.5:400,[r0 0 0 v0(i)/r0]);
    [x,y]=pol2cart(u(:,3),u(:,1));
    X= [flipud(x);x];
    Y=[-flipud(y);y];
    % text(-100-(i-1)*100,0-(i-1)*300,tel(i));
    hold on
    plot(X,Y);
    hold off
end
xlabel('x');
ylabel('y');
legend('m','E < 0', 'E = 0', 'E > 0');
%%
T=150;k=5000;E=-25;  
r0=[10 20 30 50 100 200 250];     
v0=sqrt(2*E+2*k./r0);  
F=@(t,u)[u(2);u(1).*u(4).^2-k.*u(1).^-2;u(4);-2.*u(2).*u(4)./u(1)];
plot(0,0,'r.');
hold on
for i=1:length(r0)
    [t,y]=ode45(F,0:0.01:T,[r0(i) 0 0 v0(i)/r0(i)]);
    [X,Y]=pol2cart(y(:,3),y(:,1));
    plot(X,Y);
end
function ydot=fffunction(t ,y)
    global G m
    ydot= [y(2);y(1)*y(4).^2-G*m/(y(1)^2);y(4);-2*y(2)*y(4)/y(1)];
end

